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Abstract 

We use a stochastic metapopulation model to study the com- 
bined effects of seasonality and spatial heterogeneity on dis- 
ease persistence. We find a pronounced effect of enhanced 
persistence associated with strong heterogeneity, intermedi- 
ate coupling strength and moderate seasonal forcing. Ana- 
lytic calculations show that this effect is not related with the 
phase lag between epidemic bursts in different patches, but 
rather with the linear stability properties of the attractor that 
describes the steady state of the system in the large popula- 
tion limit. 



1 Introduction 

Demographic stochasticity due to the probabilistic nature of 
events such as births, deaths, mating and disease transmis- 
sion, plays a major role in the dynamics of small popula- 
tions. Its impact was acknowledged more than fifty years 
ago by Bartlett [2|, who introduced the concept of critical 
community size. Originally defined as the population size 
above which the expected time to fade out after an epidemic 
exceeds a certain period, it is usually taken in more gen- 
eral terms as the threshold population for (a given definition 
of) disease persistence. It became a central concept in epi- 
demiology, much revisited in several attempts to provide less 
arbitrary definitions and to reconcile theoretical estimates 
with data fl8l l24l . Threshold levels of host abundance are 
equally important in ecology, a context in which the idea of 
the stochastic Allee effect was introduced 11201 to represent 
demographic stochasticity. 

The fact that many natural populations experience annual 
abundance troughs establishes an obvious connection be- 
tween average population size and extinction probability, on 
one hand, and seasonality, on the other (32). Indeed, the 
annual and multiannual incidence patterns of many infec- 
tious diseases show that seasonality is a key ingredient in 
the overall dynamics of these diseases. Despite the mathe- 
matical difficulties involved, theoretical studies have there- 
fore tried to take seasonality into account ever since the ear- 
liest efforts 1128) . The complex interplay between seasonal 



forcing and the system's nonlinearities is nowadays reason- 
ably well understood, setting the stage for the additional 
layer of complexity that arises from demographic stochas- 
ticity GHHEa. 

Another key ingredient for population persistence is spa- 
tial structure and heterogeneity. Spatial structure was first 
addressed using reaction-diffusion equations that success- 
fully modelled the spread of the epizootic in animal borne 
diseases (16) • In these models, inspired by physical systems, 
the interactions are local and the population is distributed on 
a plane. More recently, developments that explore the role 
of individual mobility and long range interactions have come 
up in the form of metapopulation models, where a number of 
typically weakly interacting units represent well mixed ho- 
mogeneous population patches JTTl[T3ll2Tll22ll25l [Tl. Along 
standing idea associated with the concept of metapopulation 
is that persistence is favoured in a fragmented population, 
provided that movement between patches accompanies spa- 
tial dispersion (4] [14). This idea has recently been shown to 
be less straightforward than previously thought lfi"2"][l5l . 

Among many aspects treated in these studies on spatially 
extended systems, the degree of synchrony of population 
abundance oscillations has received special attention as it 
has been considered the main determinant of persistence 
JlOl l3l [8l . With few exceptions associated with chaotic os- 
cillations, it has been found that a small amount of coupling 
between population patches is enough to induce synchrony 
(6] EH |9l Although they are usually called spatially 
heterogeneous, the metapopulation models in these studies 
assume the same, or very similar, parameter values for the 
different population patches, and we will refer to them as 
'uniform', keeping the term 'heterogeneous' for extended 
systems that include significant parameter variation across 
different patches. In line with available data for large urban 
populations OTl l23l . the synchronized oscillations in uni- 
form systems are moreover found to be in phase between 
patches, or, in the case of the 2-year cycle typical of measles, 
in phase opposition. This is in contrast with the results for 
heterogeneous systems, where synchronous states may cor- 
respond to intermediate phase lags If3~l l27l . 

In this paper we present an extensive computational study 
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of the combined influence of seasonality and heterogeneity 
on disease persistence. The basic unit of our model, which 
we will call a 'city', is formed by a number of individuals 
undergoing well-mixed stochastic infection dynamics whose 
parameters are specific to that city and may present seasonal 
variation. The number of individuals that interact in this way 
may comprise commuters from another city, as well as the 
residents of that city. Disease persistence is measured over 
sets of stochastic simulations of the model. We find that it 
depends in a nonintuitive way both on the level of season- 
ality and on the magnitude of the flow of commuters. We 
also find that the epidemic phase lags generated by city het- 
erogeneity have no significant effect on disease persistence. 
An analytic description of the incidence fluctuations based 
on van Kampen's expansion shows that this increase in per- 
sistence is instead related with the stability properties of the 
attractor that describes the steady state of the system in the 
large population limit. 

2 Methods 
2.1 Model 

In this section, we briefly present the metapopulation 
susceptible-infectious-recovered (SIR) stochastic model in- 
troduced in ||26l |27l to describe several interacting cities, 
which are population patches where interactions between in- 
dividuals are taken to be well mixed. 

The SIR model consists of three classes of individuals: 
susceptibles, infected and recovered. We denote their num- 
ber among the residents of city k by Sk, Ik, Rk> respec- 
tively. These numbers change due to birth, death, infection 
and recovery, which in the stochastic version of the model 
are taken as stochastic events with certain rates. As usual 
when working with time scales for which there are no major 
demographic changes, we assume that the number of indi- 
viduals that reside in city k, Nk, is fixed, so that Sk and 
Ik together completely determine the state of city k. The 
birth/death rate fi is taken to be constant, and infected indi- 
viduals recover also at a constant rate 7. When a given dis- 
ease spreads in a city, the rate of infection is proportional to 
the number of encounters between susceptibles and infected 
that take place in that city, which in turn, assuming that in the 
city the population is well mixed, is proportional to the prod- 
uct of the number of susceptibles and the number of infected 
in that city. Now these numbers should take into account 
the flow of commuters from and to that city. In the simplest 
version of the model, we will assume that the coupling be- 
tween cities 1 and 2 may be described by a single parameter, 
/, which is the fraction of the number of residents of each 
class of city 1 (respectively, 2) that are present in city 2 (re- 
spectively, 1) at any given time. The parameter / must be 
interpreted as the overall fraction of time that an individual 
from one city spends in the other city, averaged over all types 
of stays with their typical frequencies and durations. In gen- 
eral, / should be taken class and city dependent (see Supple- 
mentary Material (SM)), but we will explore here only the 



simplest case. 

The usual SIR rate of infection then becomes, for suscep- 
tible residents of city 1 while in city 1, 

[(l-/)/i+// 2 ]/M l5 

where fli is a parameter that reflects the urban characteristics 
of city 1 through the rate of encounters they elicit, and Mi = 
(1 — f)Ni + fN2 is the number of individuals present in 
city 1 at any given time. The rate of infection of susceptible 
residents of city 1 while in city 2 will be given by 

MSt [(1 - f)h + fh] /M 2 , 

with A/2 = (1 — f)N2 + fN\. Similar expressions hold for 
the rates of infections taking place in city 2. 

Our mechanistic model thus leads to represent the in- 
teraction between population patches as a weighted distri- 
bution of their respective forces of infection. Along with 
other metapopulation models based on a description of the 
underlying mobility patterns |fT9l [T7l . it extends the tradi- 
tional phenomenological modelling of interacting population 
patches by means of a single coupling parameter [22], with 
the important difference that the parameters (3k are allowed 
to differ from patch to patch, so that spatial heterogeneity 
does not come from 'patchiness' of the population only. 

The parameter fik may be time dependent to represent sea- 
sonal variability of social intercourse, or of other ingredients 
such as for instance weather conditions that influence the rate 
of infectious contacts. We will consider a time dependence 
of the form (3k(t) — /3j?(l + ecos27ri), where t is the time 
measured in years and e represents the amplitude of seasonal 
forcing. More realistic forcing terms that include a represen- 
tation of school term calendars are commonly found in the 
literature on childhood infectious diseases (e.g. [18]), but 
we expect the overall picture revealed by varying (3% and e to 
be largely independent of the particular form of the periodic 
forcing. 

With these assumptions, the stochastic process is gov- 
erned by the master equation for the time evolution of P n (t), 
the probability distribution for finding the system in state n 
at time t Poll: 

(i) 

where n denotes the state of the system given by the numbers 
of infected and susceptibles in each city and T Q (n|n'), are 
the (possibly time dependent) transition rates from the state 
n' to the state n that result from the birth-death, recovery 
and infection processes. These rates are given explicitly in 
theSM. 

2.2 Theoretical Analysis and Simulations 

A deterministic description of the model leads to a set of or- 
dinary differential equations (ODEs) for the evolution of the 
fractions Sk = Sk/Nk, %u — Ik/Nk of susceptible and in- 
fected individuals in each city. The behaviour of the stochas- 
tic system approaches this description in the limit of infinite 
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population sizes where fluctuations can be neglected. Oth- 
erwise, for large but finite systems s k and i k will fluctuate 
around the solutions of the deterministic ODEs. Using van 
Kampen's system size expansion ||30l , these fluctuations are 
approximately described by Langevin equations, which for 
two cities can be written in compact form as 

dz 4 

-j£ = Y, J J KZK + r lj{t), J = l,...,4, (2) 

at K=l 

where z = (zi, z 2 , z 3 , Z4) = (sci, x 2 , yi, 2/2), x k and y k are 
the relative fluctuations of the number of susceptible and in- 
fected in city k, and rjj(t) are zero mean Gaussian noise 
terms. The derivation of these analytic approximations, as 
well as explicit expressions for its coefficients and for the 
noise correlation function, are given in the SM. 

Neither the full stochastic system, Eq. (Q}, nor the de- 
terministic ODEs can be solved analytically. The latter are 
however amenable to qualitative analysis, which we will use, 
together with numerical integration, to study their attractors 
and their stability. As to the former, Eq. d2) can be used 
to compute the state variables fluctuation power spectrum 
and phase spectrum [27 1, which determine the amplification 
(the overall power of the fluctuation spectrum), the coher- 
ence (the fraction of the total power in a small frequency 
range around the dominant frequency), and the phase lag be- 
tween cities. Details are given in the SM. 

The simulations implement the stochastic process de- 
scribed in the preceding section on the state variables 
(Sk,Ik), using Gillespie's algorithm adapted to account for 
time dependent reaction rates [7 |. We compute the average 
extinction time (AET) for a metapopulation as the average 
over 10 3 simulations of the time it takes for the system to 
reach a state with no infectives, starting from different initial 
conditions near the equilibrium of the deterministic coun- 
terpart of the model. The simulation time was taken long 
enough for all runs to eventually become extinct. 

The crucial epidemiological parameter Ro of the unforced 
deterministic SIR model, the basic reproductive ratio, is 
given by Ro = /3/("Y + A*)- It is well known that Ro = 1 
is the critical value that separates the trivial regime, where 
the disease dies out, from the endemic regime, where the 
disease persists for ever. We will take Ro, e and / as the free 
parameters in this study. The parameters \i and 7 are kept 
fixed at 1/50 yrand 1/13 d. 

3 Results 

To illustrate the nontrivial interplay between disease trans- 
missibility and amplitude of seasonal forcing we start in this 
section by studying the stochastic SIR model in a single city. 
The results are shown in Figure 1. For e — 0, the AET 
increases exponentially with Rq, but a more complex depen- 
dence on Ro is found when e > 0. More specifically, for 
e = 0.05 the AET gets larger as we increase Ro from 12 
to 17 where it attains a maximum. A further increase of Ro 
reduces the AET, which reaches a minimum at Ro rs 21 and 
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Figure 2: Deterministic biennial attractor for the 1-city SIR 
model. 



starts to increase again from then on. For stronger season- 
ality, e = 0.12, the curve describing the dependence of the 
AET on Ro has only one maximum at about Ro = 15. A 
similar effect was found in [5 1 for the dependence of persis- 
tence on birth rate. 

As in 11231 . a simple qualitative explanation of the behav- 
ior seen in Figure 1 can be given in terms of attractors of 
the deterministic model and their stability properties. For 
the unforced case, the attractor is a stable fixed point with 
nonzero densities of susceptible and infective individuals. 
As Ro increases, the infective density associated with this 
attractor increases, as well as its stability. As a consequence, 
in the stochastic model the relative amplitude of the oscilla- 
tions of the number of infected around their average value 
gets smaller, causing the AET to increase with Ro- In the 
presence of seasonality, e > 0, the attractors of the sys- 
tem are limit cycles with periods multiples of 1 year. De- 
pending on the values of Ro, annual or biennial cycles are 
observed. Specifically, for e = 0.12 and R G [12,14.5] 
or Ro £ [15, 24], the cycles are annual or biennial, respec- 
tively. The period doubling bifurcation occurs in the interval 
14.5 < Ro < 15 that separates the regions of increasing and 
decreasing extinction times. As shown in Figure 2, the de- 
terministic biennial cycle changes with increasing Ro so that 
the infective density stays low for a longer fraction of the pe- 
riod. Due to this, in the biennial regime the AET decreases 
with Rq. In the annual regime, the stability of the cycle and 
the infective density averaged over the cycle increase with 
Ro, explaining the increasing AET for Rq € [12, 14.5] as in 
the unforced case. A similar analysis holds for e = 0.05. 

We now explore the situation where different population 
patches represent very different urban environments, as may 
be the case of an active city centre and its suburban satellite 
towns. We will therefore consider our model for two cities 
with different values of Ro, taking for simplicity equal pop- 
ulation sizes and the symmetric coupling described in Sec- 
tion l2.ll The two values i?oi an d R02 that characterize SIR 
dynamics in each patch are given in terms of a single pa- 
rameter 6 as 18 ± S, thus ensuring that the average disease 
infectiousness is kept constant as the spatial heterogeneity 
increases with 8. The qualitative results are independent of 
this assumption (results not shown), which however reduces 
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Figure 1: AET from simulations for the 1 -city SIR model with N = 4 x 10 5 (results for other population sizes are given 
in the SM, Figure SI). The arrows highlight the parameter regions where the attractor of the forced system changes from 
annual to biennial (black and red arrows), and from biennial to annual (grey arrow). 
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Figure 3: Results for the 2-city SIR model with Ni = N2 = 2 x 10 5 . (a)-(c) AET from simulations for different levels of 
seasonality, e. (d)-(f) Analytical results for the amplification and coherence in one of the cities (similar results for the other 
city are shown in the SM, Figure S2) and phase lag between cities for e = 0. 
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Figure 4: Typical time series of new disease cases per week for the 2-city SIR model with Nx = N2 = 2 x 10 5 and e = 0. 
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Figure 5: AET from simulations for the 4-city SIR model in which a central city with Ni = 2.1 X 10 5 is surrounded by 3 
non-interacting satellite cities with N2 = N3 = N4 = Ni/3. Fraction of commuters from each of the suburban cities to the 
central city, /, is 10 times larger than that from the central city to that suburban city. Analytical results for the amplification, 
coherence and phase lag of the fluctuations are given in the SM, Figure S3. 



the number of free parameters and is consistent with the fact 
that estimates of Rq come from coarse grained data. We will 
study the three levels of seasonality e — 0, 0.05, 0.12 con- 
sidered in the study of one city, and the free parameters for 
each value of e are the coupling strength / and 8. 

The results found for the AET are shown in Figure 3, top 
row. For the highest level of forcing, e = 0.12, AETs are 
very short and the dependence on / and S is weak, with only 
a slight increase in persistence for small coupling strengths, 
independently of heterogeneity. This is in line with results 
reported in the literature for uniform systems [fl4l . and can 
be understood as temporal heterogeneity superseding spatial 
heterogeneity for this level of seasonality, bringing the sys- 
tem close to extinction for a large part of the year. 

The picture that emerges for e = and e = 0.05 is signifi- 
cantly different. Extinction events are rare, and when spatial 
heterogeneity, measured by 5, is small, the AET increases 
with the coupling strength /. This is to be expected, be- 
cause disease that goes extinct in one of the identical patches 
is more likely to be reignited by surviving infectives in the 
other patch as / increases. More surprisingly, a second and 
dominant effect of persistence enhancement shows up for in- 
termediate coupling strengths when the degree of spatial het- 
erogeneity is high. 

In order to understand the causes of this remarkable in- 
crease in persistence, we have analyzed the amplification, 
the coherence and the phase difference for e = and all 
the parameter values used to compute the AETs. The results 
are shown in Figure 3, bottom row. In the left and middle 
panels, it can be seen that the enhanced persistence effect 
is associated with a pronounced reduction both of the over- 
all power of the infective number fluctuations and of their 
coherence. In the right panel, we see that the phase lag 
between patches introduced by spatial heterogeneity, which 
would seem a good candidate to explain the observed effect, 
has no significant bearing on persistence. 

In Figure 4 we plot two typical time series taken for 6 = 6 
and values of / that correspond to extreme values of amplifi- 



cation and coherence in Figure 3. These plots illustrate how 
the amplification and coherence measures translate into the 
amplitude and structure of the fluctuations of the infective 
time series, and therefore into persistence as well. Using the 
theory developed in [27 1, we show in the SM how this effect 
is associated with the dependence on / and S of the real part 
of one of the eigenvalue pairs of the linear approximation of 
the deterministic system at the endemic equilibrium. 

Finally, in Figure 5 we show a plot of AETs similar to 
that of Figure 3, but in this case for a configuration of four 
population patches, one central city connected to three non- 
interacting satellite suburban areas with equal population 
sizes. The results shown correspond to the central city and 
to one of the suburban areas for the case when the fraction 
of commuters from each of the suburban cities to the central 
city, /, is 10 times larger than that from the central city to 
that suburban city. Both the phenomenological description 
and the theoretical interpretation given above for two cities 
hold in this case as well. 

4 Discussion 

Using a stochastic SIR metapopulation model with one pop- 
ulation patch connected to one or several non-interacting 
patches, we have explored the combined effects on disease 
persistence of seasonality and spatial heterogeneity. The 
parameters in this analysis are the level of seasonality, the 
strength of the coupling, and the degree of heterogeneity 
measured by the difference between patches in the rates of 
potentially infectious contacts. We have considered only the 
simplest metapopulation structures and couplings. Results 
in more general settings (not shown) indicate that the effects 
described in this paper are robust with respect to different 
choices of these interaction parameters, and also to changes 
in the average Rq. 

In contrast with spatially structured systems formed by 
similar patches, the inherent heterogeneity of the model in- 
duces well defined phase lags in the epidemic bursts that take 
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place in different patches, suggesting a simple mechanism 
through which heterogeneity might contribute to an increase 
in disease persistence. However, we have found no clear ev- 
idence of such relation. This negative result can be under- 
stood because epidemic bursts come in short spikes, after 
which the system remains for a relatively long time close to 
extinction. The overall duration of the regime characterized 
by very low population numbers in two interacting patches 
is only slightly reduced by the phase lag. We speculate that 
in systems where the population fluctuations are more regu- 
lar a relation between persistence and the phase lags due to 
spatial heterogeneity would be found. 

We find instead a remarkable effect of enhanced persis- 
tence associated with strong spatial heterogeneity, interme- 
diate coupling strength and moderate seasonal forcing. We 
should point out that the figures illustrating this effect are 
shown in logarithmic scale for the coupling strength /, and 
that therefore the range of values of / that produce enhanced 
persistence is quite small. For these values however the ex- 
tinction times raise very significantly. Using well known the- 
oretical methods, this effect can be related with an increase 
in the stability of the attractor that describes the system for 
very large population sizes. 
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Supplementary Material 



Impact of commuting on disease persistence in heterogeneous 

met ap opulat ions 

G. Rozhnova, A. Nunes and A. J. McKane 

In this supplementary material we will specify the model we are studying and outline 
the methods we use to study it. In the final section we will also give some additional 
{Nj ' results. Our aim in the first three sections is to give a non-technical introduction to the 
methods we use; technical descriptions already exist in several of our published papers 
[1-4]. In addition, [5] gives a non-technical review of stochastic modelling and [6] a more 
mathematical treatment of the same topic. We will explain the formalism in detail for 
l/-) . the usual (one-city) model in such a way that the generalisation to two, or n, cities is 
CN \ straightforward. 

S ■ 1 The Stochastic SIR model 



> 



1.1 Formulation of the individual-based model and the master 



O 

^-i equation 

We begin with the simple SIR model — or in the context of this paper, the one-city SIR 
model. We define it in terms of "reactions" between the N constituents of the system, 
which in this case are N individuals which are either of type 5* (susceptible), type / 
(infected) or type R (recovered). These individuals are born and die at the same rate, 
and in addition these events are linked, so that the total population N remains constant. 
Thus instead of having distinct birth and death events, only combined birth/death events 



occur. Since all newly-born individuals are susceptible, these events correspond to an 
infected individual being replaced by a susceptible individual or a recovered individual 
being replaced by a susceptible individual. A susceptible individual dying results in no 
change in the state of the system, since they are simply replaced by another susceptible 
individual. 

rS ■ The reactions that define the one-city SIR model are then: 

1. Infection S + I -A I + I 

2. Recovery / — > R 

3. Birth/death I -A S 

4. Birth/death R -A S. 

Here /3, 7 and \i are respectively the rate of infection, the rate of recovery and the rate of 
birth/death. The reactions will be labelled by a (= 1,2,3,4). 

Frequently, birth and death processes are assumed to happen at the same rate, but 
remain as distinct events. This still results in fluctuations in the total population size, 
even if the average population size remains constant. By linking these events at the 
stochastic level, the population size remains constant at any system size, so that we can 
still eliminate the variable relating to recovered individuals: R = N — I — S . This means 
that there are only two variables which define the state of the system: the number of 
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susceptible individuals at a given time which we will denote by S and the number of 
infected individuals at a given time which we will denote by /. The use of the same 
symbol for an individual of a particular type and the number of individuals of that type 
should cause no confusion. The general state will then be denoted by n = (S, I). 

The four processes (JT]) in the SIR model which cause transitions to a new state have 
the following transition rates: 

1. Infection T^S - 1, 1 + 1\S, I) = 0SI/N 

2. Recovery T 2 (S, I — 1\S, I) = 7/ 

3. Birth/death T 3 {S + 1, 1 - 1\S, I) = fjtl ^' 

4. Birth/death T A (S + 1, I\S, I) = /i(N-S-I). 

We use the convention whereby the initial state is on the right and the final state is on 
the left, so that T(n|n') represents the transition rate from n' to n. 

The transition rates (T5]) define the individual-based SIR model, and when substituted 
in the master equation (Eq. (1) of the main text), give an equation for the probability, 
P n (t), of finding the system in the state n at time t. Since the master equation can only 
be solved for very simple linear systems, we require an approximation scheme to make 
progress. Before we discuss this, we first derive the macroscopic description of the SIR 
model (in the form of differential equations) from the master equation. 

To do this it is useful to write the master equation in a slightly different form, by 
introducing the stoichiometric coefficients u a , for reaction a. They indicate by how much 
S or / increase or decrease in a given reaction, that is, n. = n' + u a . So for instance in 
reaction 1, the components of V\ are = — l,u^ = +1, and for reaction 2 we have 

v% = 0)^2 = ~ 1) an d so on. Here the superscript (1) refers to S and (2) to /. Then 
the master equation can be written as 

dP (t) - 

-^f 1 = [ T «( n l n - v«)Pn-v a (t) - T a {n + u a \n)P n (t)} , (3) 
where in this case M — 4. 



1.2 The macroscopic equation 

The macroscopic equation is found by multiplying Eq. ([3]) by n, and summing over all 
possible values of n. After making the change of variable n — > n + u a in the first 
summation, one finds that 

^P = E"«< r «( n +^»>> ( 4 ) 

a=l 

where the angle brackets define the expectation value: 

(...} = J2(---)P n (t). (5) 

n 

The first approximation we will make (which is exact in the limit N — > 00) is to 
label the states of the system by the fraction of the population which are susceptible and 
infected, rather than the numbers of susceptible and infected individuals. In other words, 
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we will use s = S/N and i = I/N (in the limit iV — > oo) rather than S and /, and treat 
these as continuous variables. This is the sometimes called the diffusion approximation. 
So the state of the system will be determined by the new variables 

d> = (s,i) = lim - — . 
Dividing Eq. (j3j) by N, we have that 

d<fi M 

-fa =^2 U a9a(<j>), (6) 

o=l 

where the function g a {4>) is defined by lini7v->oo N^ 1 (T a (n + v a \n)). 

Equation (jfJJ) is the usual deterministic (macroscopic) equation for the SIR model, and 
specifies how the mean fraction of susceptible and infected individuals changes with time. 
To write it in a more familiar form, we first note that in the limit iV — > 00 the system 
is deterministic, and so the variance (and higher cumulants) of n are zero. This implies 
that 

lim N~ l (TJn + v a \n)) = lim N' 1 T a ((n) + v a \(n)). 
The mean value (n) can now simply be replaced by N<f>, leading to the identification: 

g a {4>)= lim N- x T a {N$ + v a \N<l>). (7) 

N— >oo 

For example, from Eq. (J2J, T X (S - 1, / + 1\S, I) = 0SI/N, and so g l = lim^oo A^ 1 7i = 
[i si. This leads to the following characterisation of the reactions — which is the only 
information we require for specifying both the deterministic and stochastic dynamics of 
the system: 

1. Infection 

2. Recovery 

3. Birth/death 

4. Birth/death 

We define 

M 

A(<J>) = J2»a9a(<j>), (9) 

a=l 

so that the macroscopic equation (JB]) is 

g-A(„. {10) 

Writing this in terms of the components of <fi = (s, i) we find, using Eq. (JBJ), the well-known 
equations for the deterministic SIR model: 

— = — psi + 111 + a ( 1 — s — i) , 
at 

di 

— = +(3si — 72 — iii. (11) 





= -1, 


(2) 


+ 1 


9i 


= (3si 




= 0, 


u 2 — 


-1 


92 


= 7« 


"3 


= +1, 


(2) 


-1 


93 


= fii 


"4, 


= +1, 


(2) 





94 


= /i(l 



3 



As is also well-known these equations have two fixed points: one where there are no 
infected individuals in the population, that is, (s*,i*) = (1,0), and the one which is of 
most interest: 

* 7 + A* „•* A* \P ~ (7 + »)] 



' Pil + fi) 

In the above the asterisk denotes a fixed point of the differential equations. 



(12) 



1.3 The stochastic equation 

A question of immediate interest which can be answered from the deterministic equation 
( ITT1) is the long term behaviour of the model, which involves finding the stability of the 
fixed points. We carry out this analysis in the usual way, by linearising the equations ( ITT]) 
about the fixed point. That is, we substitute 

s = s* + 4=, i = i* + -^=, (13) 



into the equations fill I) , keeping only linear terms in either x or y. This leads to the 
equations 

dx 

-T- = Jnx + J 12 y, 
at 

^| = Ji2x + J 22 y, (14) 

where the Jacobian J is evaluated at the non-trivial fixed point 

/-P?-li -/3s* \ ( "(7 + ^)\ 

V 0i* /3 S *-( 7 + /i)/ \I££+M / 

The eigenvalues of this matrix have a negative real part, and are in fact a complex 
conjugate pair for all realistic values of the rates /3,7 and \i. Thus this fixed point is 
stable and small perturbations undergo damped oscillations towards (s*,i*). 

The stochastic version of the model can now be easily written down. We will only 
be interested in fluctuations about the stationary state of system, that is, about the 
fixed point ( |T2l . Furthermore we will only investigate linear fluctuations, since numerous 
studies have shown that this captures the actual dynamics very well indeed. This is 
known as the linear noise approximation (LNA) or the next-to-leading order in the van 
Kampen system-size expansion. It consists of implementing the expansion ( TT3|) but on 
the master equation, not on the deterministic equations. The latter calculation, described 
above, amounted to an analysis of what happened to a single small perturbation to the 
deterministic dynamics. The linearisation carried out on the master equation, in contrast, 
results in extra "noise" terms which are added to (114j) and which continuously modify the 
(stochastic) perturbation. In this context the factors of iV _1//2 in the linearisation (|T3|) are 
important; they encode the fact that fluctuations will naively be of this order, and so pick 
out the linear correction in the expansion procedure. The N~ l l 2 factors were included 
in the linear stability analysis of the deterministic equation to highlight the connection 
with the LNA, but clearly have no effect on Eq. (fT4l . since any factor multiplying the 
perturbation will cancel out in this linear equation. 
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We will not carry out the calculations leading to the equations describing the stochastic 
fluctuations about the fixed point, since they have already appeared several times in the 
literature in various forms. The equations are 

dx 

— = Jux+ Jizy + mit), 
at 

dy 

= Jl2X + J 22 y + V2(t) 7 (16) 

where r]i(t) and r] 2 (t) are noises with a Gaussian distribution and zero mean, and with a 
correlation function given by 

(r)j(t)vK(t'))=B JK 5(t-t'). (17) 

Here Bjk is the matrix 

M 



B JK (4>) = ^ »i K) 9a(4>), J,K = 1, 2, 



a=l 



evaluated at the fixed point f|T2|) . Using Eqs. (jSJ) and ffT2l one finds that 

/flrt-. + rti-O -0s;- x / SdtteM - gteMJfctell \ 
V -Bs*i*-ai* Bs*i* + (j + u)i* J \ _Mr±MMzto)l 2 M [/3-( 7+M )] 

(19) 

Equations (|T6l) . (JTTJ) and (|T9l completely specify the stochastic fluctuations about the 
deterministic SIR model. 

To analyse the stochastic differential equations ( JT6l) further it is convenient to take 
their Fourier transforms and introduce 

/OO /*00 
e iu}t x(t)dt, y(u)= e iult y(t)dt. (20) 
-oo J —oo 

The lower limit of the integration comes from the fact that to ensure that the system is 
in a stationary state, the initial conditions have to be set in the infinitely distant past. 
Since the Fourier transform of dx/dt is —iux, one finds that 

[-iui - Jii] x(u) - J\2y(u) = ^i(w), 
-J2ix{u) + [-zw - J22] y{u) = fj 2 {co). (21) 

Defining the matrix —iojSjk — Jjk to be <&jk{w), Eq. ( 1211) may be written in the more 
concise form Y^k=i &jk{u)z k (uj) = fjj(u), where z\ = x and z 2 = y. This may be solved 
to yield 

2 

~zj(u) = J2<S>jk(u)Vk(u), (22) 

K=l 

where $ _1 is the inverse of the matrix $. 
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The nature of the stochastic fluctuations can usefully be investigated by calculating 
the power spectrum defined by 

2 2 

PM = (I^HI 2 ) = J2J2®JK^)Bkl (23) 

K=l L=l 

where <&t is the Hermitian conjugate of the matrix $. Examples of power spectra from 
a number of models are given in [1-4]. To characterise different spectra we use two 
measures. The amplification is defined as the overall power of the fluctuation spectrum. 
The coherence is the fraction of the total power in a 10% frequency range around the 
dominant frequency of the spectrum. 

The explicit form of the denominator in Eq. (1231) is |det$(u;)| 2 . From this we can 
deduce that a pair of complex conjugate eigenvalues of J with small real parts will give 
rise to a peak in the power spectrum, whose location in frequency is determined by 
the imaginary parts [3]. The file ESM.mov included as Supplementary Material shows 
how the persistence enhancement effect discussed in the main text is associated with the 
dependence of one such pair of eigenvalues of J on the parameters of the system. 

2 The Two-City Stochastic SIR model 

Having discussed in detail the formalism and general procedure in the one-city case in a 
way that naturally generalises, the description of the two-city case can be more concise. 
What is novel is the nature of the interactions between the two cities, which is what we 
discuss first of all. 

2.1 Form of the interaction between the two cities 

The interaction between the two cities is only reflected in the first type of reaction in 
Eq. ([I]); the other three remain essentially unchanged since they involve only one individ- 
ual, and so only one city. 

The number of individuals in the three classes belonging to city j (j = 1,2) will be 
denoted by Sj , Ij and Rj respectively. We again assume that births and deaths are coupled 
at the individual level, so that when an individual dies another (susceptible) individual 
is born. Therefore the number of recovered individuals is not an independent variable: 
Rj = Nj — Sj — Ij, where j = 1,2. Since we do not focus on specific individuals, we 
will not be concerned with the precise movements of commuters between the cities - 
the frequency or duration of their commute — but only the fraction of the population 
which is away from its home city at any given time. For commuters from city j to city k 
we will denote this fraction by It follows that the number of individuals in city 1 is 
M 1 = (1 - / 21 )Aq + f 12 N 2 and in city 2 is M 2 = (1 - f 12 )N 2 + f 21 N x . 

We will assume that the birth/death rate and the recovery rate are the same in both 
cities, but that the infection rates are city dependent: (3i in city 1 and f5 2 in city 2. When 
including seasonal forcing, as we are in this study, the infection rates are time-dependent: 
f3j(t) = fl® (1 + e cos 27r£), where the time t is measured in years and e represents the 
amplitude of the seasonal forcing. 

There are four different types of infection events. 

(i) Infective residents of one city (say j) can infect susceptible residents of the same 
city. The rate for this to occur is (3j (1 — fkj) Sj (1 — /jy) Ij/Mj. 
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(ii) Infective commuters can travel from city k to city j and infect susceptible residents 
in their home city j. The rate for this to occur is f3j (1 — /jy) Sj fj^Ik/Mj. 

(hi) Infective residents in their home city (say k) can infect susceptible commuters from 
the other city (j). The rate for this to occur is (3k fkjSj (1 — fjk) Ik/M k . 

(iv) Infective commuters can infect susceptible commuters away from their home city of 
j. The rate for this to occur is (3k fkjSj f k jIj/M k . 

Adding these rates together we obtain the total transition rate for infection of Sj 
individuals as 

+ Pik^, (24) 



N, 



where 



,/, ; ) 2 .V f fafijNj _ ft (1 - f kj ) f jk N k 4,/, f il ,/- / j.V/, , . 



Ma 



Equation (|24|) gives us the generalisation of the infection rate to the two-city case. 
There are now eight processes rather than the four processes ([1]) of the one-city SIR 
model which cause transitions from one state to another. They are 



1. 


Infection of Si 






S,. 


h) 


a S-Lh 

- Pll— 


+ /3l2 


Sih 

N 2 


2. 


Infection of S2 




T 2 (S 2 -1,I 2 + 1 


s 2 . 


h) 


_ a S2I1 

- P211vT 


+ P22 


S2I2 
N 2 


3. 


Recovery of l\ 




T 3 (Sx,h-l 


s h 


h) 


= ih 






4. 


Recovery of I 2 




T 4 (S 2 ,I 2 -1 


s 2 . 


h) 


= ih 






5. 


Birth/death in 


city 1 


T 5 (3! + 




h) 


= fill 






6. 


Birth/death in 


city 2 


T 6 (S 2 + 1,I 2 -1 


s 2 . 


h) 


= 






7. 


Birth/death in 


city 1 




s 1 . 


h) 


= li{N x - 


Si- 


h) 


8. 


Birth/death in 


city 2 


T 8 (S 2 + 1,I 2 


s 2 . 


h) 


= 


s 2 - 


h). 



(26) 



Note that we have not listed all the state variables as arguments of the transition rates 
T a {- • • I • • • ) — only those which are most relevant to the reaction under consideration. 

We may now set up the master equation in the same way as for the one-city case. 
Introducing n = (Si, S 2 , Ji, J 2 ), the master equation takes the form where now M — 8, 
and the reactions with the transition rates are given by Eq. ( )26l) . The variables relevant 
for the deterministic and stochastic equations are Sj = Sj/Nj and ij = Ij/Nj (j = 1,2), 
with <fi = (si, s 2 , i\, i 2 ). The analogue of Eq. (JSJ) for two cities (although now with only 
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the non-zero stoichiometric coefficients listed) is 



Infection of Si 
Infection of S 2 
Recovery of Ii 
Recovery of I 2 
Birth/death in city 1 
Birth/death in city 2 
Birth/death in city 1 
Birth/death in city 2 



v 



(i) 
1 

,( 2 ) 



,(2) 



-1 
-1 



,(3) 



5 "1 



V. 



+ 1, I/, 

+1, v, 



,( 4 ) 
2 

(3) 
3 

(4) 
4 

(3) 
5 

(1) 
6 

4" 

,(2) 



+1 


9i 


= PuSih 4- 


fil2S\l2 


+1 


92 


= P2182H + 


/322 s 2^2 


-1 


93 


= 7»i 




-1 


94 


= 7^2 




-1 


95 


= fik 




-1 


96 


= 




+1 


97 


= A*(l - si 




+1 


9s 


= M 1 - s 2 


-ia)- 



(27) 



2.2 The macroscopic and stochastic equations 

The information given in Eq. ( 1271) is again sufficient to completely specify both the macro- 
scopic and stochastic equations. 

The macroscopic equation ( TTUj) . where A(</») is defined by Eq. fl9]), when written down 
in terms of the components of <fi = (si, s 2 , U, iz), is 



dsi 

~dt 

ds 2 

~dt 

dii 

~dt 

di 2 



-PuSiii - (3 12 sii 2 + nix + (jl (1 - si - i x ) 
-/3 2 is 2 ii - f3 22 s 2 i 2 + fxi 2 + s 2 - i 2 ) 

+/3nsiii + f3 u s 1 i 2 - jix - fii h 
+(3 21 s 2 i 1 + f3 22 s 2 i 2 - ji 2 - \i% 2 . 



(28) 



It is known that for two cities (and also for n cities) that a unique non-trivial fixed 
point exists which is globally stable. The Jacobian evaluated at this fixed point is 



J 



( -Pnii ~ 012% ~ A* -P11S1 

-021*1 - P22I2 - H -021*2 

Pni*i + Pi2i*2 p llS *-( 7 + n) 

fail + 022% P21S* 2 



~Pl2S*i \ 
—022S 2 
Pl2S\ 

$22**2 - (7 + y) J 



(29) 

The stochastic fluctuations about the fixed point are obtained through a generalisation 
ofEq. 03}: 



Sj = Sj + 



l j = i ) + ^7=i J = 1,2, 



(30) 
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and keeping only linear terms in either x or y. This leads to the equations 
dxi 



dt 

dx 2 

~df 

dt 

dt 



Ji\X\ + J12X2 + J13Z/1 + J142/2 + »7i (*) ) 

^21^1 + ^32^2 + JtsVi + J2M2 + T]2{t), 

J31X1 + J32X2 + JssVi + J34Z/2 + r) 3 {t), 
JaXx + J 42 x 2 + JazVi + JmV2 + rji{t), 



(31) 



which may be written in a more compact form by introducing the vector of fluctuations 

z = (x 1 ,x 2 ,y 1 ,y 2 ): 

A 

dzj 
dt 



^JjKZK+TH(t), J=l,...,4. 



(32) 



K=\ 



Here the rjj(t) are noises with a Gaussian distribution and zero mean, and with a corre- 
lation function given by Eq. ( TT7|) . 

The B matrix f ll8p for the two-city case, evaluated at a fixed point, is 



B 



5(1) 


B m 1 


£>(3) 


5(4) 



(33) 



5(1) 



^(2) = S (3) 



where these submatrices are given by 

Pxis\i\ + Pi2sl% + //(l - si) 

ftl^ + P22S* 2 i* 2 + ~ S* 2 ) 

-Pns\i\ - (3i2s\i* 2 - ni\ 

-/3 21 s* 2 il - l3 2 2S* 2 i* 2 - fii* 2 

/9ns*i* + Aasjia + (7 + ^)*J 

ftuS^ + /9 22 S2*2 + (7 + A*)*2 

Using Eq. (I28p at the fixed point these may be simplified to 

2//(l-sJ) 



and 



5(4) 



(34) 



(35) 



(36) 



and 



5(1) 



5(2) = B (3) 



5(4) 



2/i(l - s*) 
-(7 + 2/iK 

-(7 + 2^* 

2(7 + /iK 

2(j + n)i* 2 
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(37) 



(3? 



(39) 



We now introduce the matrix 



JK 



(zj(u)z* k (uj)) 



EE 

L=l M=l 



(*t) 



MK 



(40) 



where here (and only here) * denotes complex conjugation. In the one-city case, where 
the focus is on finding the frequencies and amplitudes of the stochastic oscillations, only 
the power spectrum (when J = K) is usually analysed. When studying the model with 
two cities, we will also be interested in the cross-correlations between infection in two 
different cities, and so will also wish to calculate the cross-spectrum (when J ^ K). It is 
frequently convenient to normalise this by the relevant power-spectrum, and instead work 
with the complex coherence function (CCF) defined by 



Cj K {u) = 



Pjk(co) 
'Pjj(uj)P kk \ 



UJ) 



(41) 



The CCF will in general be complex for J ^ K, and so typically one calculates its 
magnitude and phase. The phase is given by 



u) = tan 



hn(C JK {u>)) 
Re(C JK (u)) 



tan 



lm{P JK (u)) 
Re(P JK (u)) 



(42) 



The phase used in the main text is found by evaluating (142!) at the value of u that 
maximises the modulus of the CCF (1411) [4]. 



3 The n-City Stochastic SIR model 

We will be relatively brief in this section, and only outline the results, since the formalism 
and general procedure is as for the two-city case. The main difference is that there is 
a fifth type of infection event — in addition to those mentioned for two cities. This is 
due to the fact that infective individuals can commute from city k and infect susceptible 
individuals from city j in city £, where j, k and I are all different. This is only possible 
when there are three or more cities. 



3.1 Form of the interaction between the n cities 



The number of individuals in the three classes belonging to city j are denoted by Sj, Ij 
and Rj = Nj — Sj — Ij as before, where now j = 1, . . . ,n. We will also introduce the 
notation 

so that the number of individuals in city j may be written as 



Mi 



k+j 



[l-f 3 )N 3 + Y,fjkN k . 

k+j 



(44) 
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As for the two city case, we will assume that the birth/death rate and the recovery rate 
are the same in both cities, but that the infection rate for city j is ftj. 
As mentioned above, there are five different types of infection events: 

Infective residents of one city (say j) can infect susceptible residents of the same 
city. The rate for this to occur is (3j (1 — fj) Sj (1 — fj) Ij/Mj. 

Infective commuters can travel from city k to city j and infect susceptible residents 
in their home city j. The rate for this to occur is (3j (1 — fj) Sj fj^I^/Mj. 

Infective residents in their home city (say k) can infect susceptible commuters from 
the other city (j). The rate for this to occur is j3 k fkjSj (1 — fk) h/Mk. 

Infective commuters from city j infect susceptible commuters from j in city i (£ ^ j). 
The rate for this to occur is f3i fySj fgjlj/Mg. 



m 



iv 



(v) Infective commuters from city k can infect susceptible commuters from city j in city 
£ (£ 7^ j, k). The rate for this to occur is (3e fijSj fikh/Mi. 

Adding these rates together we obtain the total transition rate for infection of Sj 
individuals as 



where 



N, 

k=i ' 



Pjj " M, ' ^ M, ' J 

3 m 

/ f ),/-/,.V;, P k f kj (l-f k )N k 

jk "' Mj M k 



+ J2^W^> j,k = l,...,n;j*k. (46) 



s, I, 



The analogue of Eq. f l26|) is 

1. Infection of Sj T^Sj - 1, Ij + l\Sj, Ij) = Y2=i Pjk% 

2. Recovery of Ij T 2 j(Sj, Ij — l\Sj, Ij) = jlj 

3. Birth/death in city j T 3 j(Sj + 1, Ij — l\Sj, Ij) = filj 



k 



(47) 



4. Birth/death in city j T4j(Sj + 1, Ij\Sj, Ij) = fi(Nj — Sj — Ij), 

although we have aggregated the reactions for convenience, so that each line represents n 
reactions with j = 1, . . . ,n. 

From this, the information required to deduce the macroscopic and stochastic equa- 
tions can be found: 



I. Infection of Sj = -1, z/j" +j) = +1 gij = YJk=i^jkSjik 

,(n+j) 



2. Recovery of Ij v\- = — 1 g 2 j = jij 

3. Birth/death in city j v^j = +1> u 3j + ^ = ~ 1 9sj = 



(4* 



4. Birth/death in city j u$ = +1 g±j = — Sj — i 
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3.2 The macroscopic and stochastic equations 

The macroscopic equation (fT0|) . where A(tj>) is defined by Eq. <^j, when expressed in 
terms of the components of <fi = (si, . . . , s n , ii, . . . , i n ) may be written down using the 
information in Eq. (I48p . One finds 



df 



fc=i 



(49) 



fc=i 



The Jacobian evaluated at the unique stable non-trivial fixed point is 

J = 

where these submatrices are given by 



" J {1) 


Jp) ' 







3, 



(i) 



J, 



(3) 



fad} + V 
Yl fort 



5, 



4 4) = ^ s *j - (7 + M) S Jk . 



(50) 



(51) 



Finally, the B matrix ( jl8l) evaluated at a fixed point has the form ( 1331) where 



5 



5 



(i) 

(2) 



6* 



jki 



B 



(3) 



^/3^* + ( 7 + /i) 



Using Eq. (149]) at the fixed point these may be simplified to 



5 



(2) 



[2/i(l - s*)] (J iJfe , 

= - [(7 + 2/i)**] 



sj? = [2( 7 + »)i*]5 jk . 



(52) 



(53) 



The cross-spectra, complex coherence function and phase spectra are as in Eqs. (I4t3]) - (l42l) . 
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4 Additional results 



In this section, we will give additional results for the 1-city (Figure 1), 2-city (Figure 2) 
and 4-city (Figure 3) SIR models. 

Figure 1 shows the average extinction time (AET) from simulations of the 1-city SIR 
model as a function of the basic reproductive ratio, R , for different population sizes, N, 
and levels of seasonality, e. The results for iV = 4 x 10 5 (middle row) were shown in 
Figure 1 of the main text. We observe that for smaller (N = 2 x 10 5 , top row)/larger 
(N = 8 x 10 5 , bottom row) populations the dependence on R Q of the AET becomes much 
less/more pronounced. 

Figure 2 shows the AET from simulations of the 2-city SIR model for different levels 
of seasonality, e, and analytical results for the amplification, coherence for both cities and 
phase lag between cities. The basic reproductive ratio is taken to be i?oi = 18 + 5 in one 
of the cities and i?02 = 18 — <5 in the the other city, where 5 = 0, 1, . . . , 6. Panels (a), (b), 
(c), (e), (g) and (h) are repeated from Figure 3 of the main text. The new panels (d) and 
(f) show that both the coherence and amplification in the first city are low in the region 
where the AET is longest, similarly to the results for the second city, shown in panels (e) 
and (g). 

Figure 3 shows the same results as Figure 1 but for the 4-city SIR model in which a 
central city with larger population is surrounded by 3 non-interacting satellite cities with 
smaller (and equal) populations. The basic reproductive ratio is taken to be i? i = 18 + 5 
in the larger city and -Ro2,o3,o4 — 18 — 5 in the smaller cities, where 5 = 0,1,..., 6. 
The results refer to the central city and one of the satellite cities. We consider the case 
when the fraction of commuters from each of the satellite cities to the central city, /, 
is 10 times larger than that from the central city to that satellite city. Panels (a), (b), 
(c) repeated from Figure 5 of the main text show that the AET is maximum at strong 
heterogeneity (large 5) and intermediate coupling strengths /. The additional panels (d)- 
(h) again suggest that this behaviour can be explained in terms of the coherence and 
amplification and that the phase lag between the cities does not play the dominant role 
in such behaviour of the AET. 
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Figure 1: AET from simulations for the 1-city SIR model. The arrows highlight the 
parameter regions where the attractor of the forced system changes from annual to biennial 
(black and red arrows), and from biennial to annual (grey arrows). 
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Figure 2: Results for the 2-city SIR model with JVi = N 2 = 2 x 10 5 . (a)-(c) AET from 
simulations for different levels of seasonality, e. Analytical results for (d)-(e) coherence, 
(f)-(g) amplification and (h) phase lag between cities for e = 0. 
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Figure 3: Results for the 4-city SIR model in which a central city with Ni = 2.1 x 10 5 is 
surrounded by 3 non- interacting satellite cities with N 2 = N 3 = N4 = Ni/3 = 7 x 10 4 . 
Fraction of commuters from each of the suburban cities to the central city, /, is 10 
times larger than that from the central city to that suburban city, (a)-(c) AET from 
simulations for different levels of seasonality, e. Analytical results for (d)-(e) coherence, 
(f)-(g) amplification and (h) phase lag for e = between the central city and one of the 
satellite cities. 
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